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Abstract 

In this work we apply a highly efficient Monte Carlo algorithm recently proposed by Newman 
and ZifF to treat percolation problems. The site and bond percolation are studied on a number of 
lattices in two and three dimensions. Quite good results for the wrapping probabilities, correlation 
length critical exponent and critical concentration are obtained for the square, simple cubic, HCP 
and hexagonal lattices by using relatively small systems. We also confirm the universal aspect of 
the wrapping probabilities regarding site and bond dilution. 
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I. INTRODUCTION 



Among the several methods for treating disordered systems and geometrical problems, 
percolation theory is certainly one of the most important. Due to its similarities with 
transitions that occur in many other systems (not only physical, but biological, social etc.), 
percolation has been used in studies within a large variety of fields. Forest fires [fj , biological 
evolution , epidemics j^, social influence p] and dilute magnetism p are only a few 
examples of the wide applicability of this theory, as well as percolation (e.g., in porous 
media) itself 0|. 

Although easily defined, percolation presents theoretical and computational difficulties. 
For instance, the percolation threshold for the site problem on a simple square lattice is not 
known exactly. Therefore, approximate solutions are necessary, and much effort has been 
dedicated in this direction. From the theoretical point of view, one can utilize mean-field 
y, ^ and renormalization group 13, [l], 12, techniques, among others. In particular, 
computer simulations constitute a powerful tool in this area, since their application to per- 
colation is simpler than for many other problems in statistical physics Typically, one 
can obtain a valid configuration by simply populating sites (or bonds) in a given lattice. To 
measure quantities of interest, such as the percolation threshold or the mean cluster size, it 
is necessary to identify all clusters in the configuration. For this purpose, many algorithms 
lave been developed, the best known perhaps being that devised by Hoshen and Kopelman 
15^ . Other algorithms, like hull-generation , can also be used, but only to answer 

some specific questions. More recently, Newman and Ziff proposed a new algorithm jisj ], 
which is general and quite efficient, both in its computational requirements and in its accu- 
racy. Although the algorithm can be used to obtain any observable of the problem, in their 
papers they have used it to calculate the so-called wrapping probabilities to investigate a 
number of aspects of the problem. For example, using exact values of the wrapping prob- 
abilities, a high-precision result for the site percolation threshold on the square lattice was 
obtained jisl . 

Unfortunately, for most lattices, one does not know the exact wrapping probabilities 
and, in using them, the problem must be tackled in a different way. One of the wrapping 
probabilities - R^i\ which is the probability that a cluster wraps around one specified axis, 
but not around the other ones - is particularly useful. In this case it is not necessary to know 



its exact value at the critical concentration for the infinite system since it has a maximum 
from which the critical point can be obtained. The method, however, is indeed capable to 
properly estimate the other wrapping probabilities which do not exhibit a maximum (in such 
cases there is just a crossing region close to the critical threshold). On the other hand, as we 
will see below, there are still some probabilities in dimensions higher than two that present, 
besides the maximum, a crossing region from where critical behavior is also achieved. 

In this work we compute the percolation threshold and the correlation length exponent, 
as well as the set of the unknown wrapping probabilities using the Newman-Ziff algorithm. 
These quantities are evaluated by employing the usual finite-size scaling as well as a cell- 
to-cell scheme jiol. After summarizing the Newman-Ziff approach in the next section we 
describe, in section III, the method which enables us to evaluate such geometrical quantities 
and we apply it to site and bond percolation on the square, simple cubic, HCP (hexagonal 
close-packed) and simple hexagonal lattices. Concluding remarks are given in the final 
section. 



II. THE NEWMAN-ZIFF ALGORITHM 



To determine the percolation transition, this algorithm uses the wrapping probability 
Rl{p)i which, for a given site (or bond) occupation p, is basically the probability that a 
cluster wraps around a system with periodic boundary conditions on a lattice of linear 
dimension L. This wrapping can, however, be defined in various manners, each with its 
own probability. For instance, on two-dimensional lattices one has: (i) and , the 
probabilities that a cluster wraps around the system in the horizontal or vertical direction, 
respectively (on a square lattice these quantities are equal); (ii) R^l\ the probability that 
the cluster wraps around the lattice either horizontally or vertically (or both); (iii) R^^ 
for wrapping in both horizontal and vertical directions; (iv) i?^^^ for wrapping around one 
specified direction but not the other one. Different lattices can allow for further geometrical 
choices for Rl{p). For example, on a simple cubic lattice, besides Rl\ we can define R^^^ 
as the probability that there exists a cluster that wraps the system in two directions, but 
not around the third one. Analogously to R^l\ we have i?^^"* for wrapping around the three 
directions. 

In order to evaluate these quantities it is necessary to generate many independent real- 



izations of the algorithm, each of them consisting of the following steps: 

(1) Initially, all sites are empty; 

(2) Sites are chosen to be occupied at random; 

(3) When a new site is added, one must check all its neighbors to verify if the new site forms 
an isolated cluster (all neighbors empty) or if it joins together two or more clusters. In the 
first case, we need do nothing. In the latter, we have to update the cluster listing. Clusters 
are stored in a tree structure, with one site of each cluster considered the root site. All sites 
in a given cluster, other than the root, have a pointer to some other site in the same cluster, 
such that by following a succession of such pointers one can ultimately reach the root. In 
order to join two clusters we simply add a pointer from the root of the smaller cluster to 
the root of the larger one; 

(4) Each time step (3) is repeated, we evaluate the quantities of interest as a function 
of the number n of occupied sites. may be any of the wrapping probabilities Rl. Let n' 
be the number of occupied sites for which percolation first occurs in a given realization. 
represents the fraction of realizations in which n' is less than or equal to n. Using all Q^'s 
so evaluated, it is possible to calculate the function Ql{p) for all possible values of p in the 
range between and 1 by a convolution with the binomial distribution 



For the bond percolation problem, we just replace sites by bonds in the above steps. 

The evaluation of the statistical errors can be done in a conventional fashion. As discussed 
in reference the standard deviation of the binomial distribution is given by 



where Ql{p) in the above equation has been taken as the mean value of the corresponding 
wrapping probability and MCS is the number of Monte Carlo steps per site. 

As an illustration, we show in Fig. ^a) and^b) the wrapping probabilities R^l^^ and R^l^ 
as a function of the concentration p of occupied sites for square lattices of different sizes. The 
exact values R^rivr) of these quantities for an infinite square system were derived by Pinson 
and Ziff Q- InRef. it has been used these exact values to obtain an estimate for 
the percolation threshold Pc- For each system size L, one determines the p value yielding a 
wrapping probability equal to exact critical value. This p value provides the estimate of p^ 




n 



(1) 




(2) 



for that L. In the critical region, one knows that the estimates converge to the threshold 
Pc according to 




FIG. 1: Wrapping probabilities i?^ as a function of occupation probability for site percolation on: 
square lattices (a) i?^^ and (b) R^^^'-, and hexagonal lattices (c) i?^^ and (d) R^^^ ■ For the square 
lattice the data have been obtained by taking 4.0 x 10^ MCS and lattice sizes L = 32, 48, 64, 96, 
128. For the hexagonal lattice the data have been obtained by taking 5.0 x 10^ MCS and lattice 
sizes L = 10, 12, 16, 20, 24. In all figures the error bars were omitted for a better visualization. 
Vertical dashed lines represent the lowest and the greatest values of p we have utilized for evaluating 
wrapping probabilities. 



p'^-Pc-L-'-'/r (3) 

For sq uare systems, using the known value 4/3 of the exponent v and ^ = 2, as obtained 
in [l^, we have p^ — Pc ~ By a finite size scaling, Newman and Ziff obtained 

Pc = 0.592 746 21(13) for the infinite system. This procedure is more complicated in higher 
dimensions. Since neither 6 nor u are known, one has to vary the scaling exponent to obtain 
a straight line. In Ref. |l8|], it was found that the estimates of the percolation threshold for 
a simple cubic lattice scale as L^^. Thus, it seems to be not so easy to determine u and 6 
separately by this method. 

One has to note (Fig. ^p) that R^l^ is different from the other probabilities, as it exhibits 
a maximum. In this case, p^ can be estimated from the position of this maximum. R^^^ 
is then of particular utility in systems for which the exact values are not known. We will 



see, moreover, that all the other wrapping probabilities can also be used to estimate the 
percolation threshold, as well as the correlation length critical exponent, on any lattice. 
There are, in addition, some quantities like in three dimensions which exhibit both a 
maximum and a crossing region. 

III. APPROACH AND RESULTS 

We have applied the Newman-Ziff algorithm to site percolation on the two-dimensional 
{2d) square, the three-dimensional simple hexagonal and simple cubic {3d) lattices, and to 
bond percolation on square {2d), simple cubic and HCP (hexagonal close-packed) lattices. 
Table HI gives the system and sample sizes used in our study. 

Before discussing the results, we analyze the standard deviation of some particular quan- 
tities. Fig. 121 shows the relative error of R^^^ for the square lattice and R^^^ for the hexagonal 
lattice as a function of the L, at the critical concentration. Apart from a strong dependence 
with small L, we clearly see that for larger lattice sizes the relative error is almost indepen- 
dent of L, as predicted by Eq. and previously stressed by Newman and Ziff [3], even 
for the three-dimensional hexagonal lattice. Similar behavior is found for other wrapping 
probabilities, other three-dimensional lattices, as well as for the bond problem in different 
lattice structures. 
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FIG. 2: (a) Relative error of -R^ , where AR^ = aji^^ for site percolation on the square lattice, 
(b) Relative error of R^^\ where AR^^^ = aji^ for site percolation on the hexagonal lattice. 



Let us now discuss the evaluation of the critical exponent, the percolation threshold 



and the wrapping probabilities. In order to get an idea of the performance of the present 
approach, we will first apply it to the problem in two dimensions where exact (or more 
accurate) results are available. From Fig. ^a) one sees that the derivative of b!"!^^ at the 
critical concentration increases as the lattice size increases. In fact, one expects that the 
maximum derivative of any wrapping probability not exhibiting a maximum scales as 



max 



L^l\ (4) 



dp 

Thus, the critical exponent v can be estimated without any consideration of the critical 
concentration pc by taking the scaling behavior of the derivatives of the thermodynamic 
quantities Rj,. They can be straightforwardly computed from relation (fT|) 



dp ^ 



(A^-n)p"(l -p) 



AT-n-l 



Ql ■ (5) 



In Fig. Elwe plot, on log^Q scales, the maximum value of the derivative of R^^^ as a function of 
system size, for site percolation on the square lattice (the hexagonal lattice will be discussed 
later); a linear fit yields u^'"^ = 1.334(4). Other quantities give independent estimates of the 
exponent (the corresponding data are too close to those of R^^^ to be distinguished on the 
scale of Fig. E}. We find z/W = 1.331(2); u^^^ = 1.339(4) and z/(^) = 1.327(1). Combining 
these four estimates we obtain u = 1.333(5), in very good agreement with the exact result 
4/3. 

Fig. m illustrates the approach for evaluating the critical concentration, as well as the 
wrapping probabilities at pc, applied again to the site percolation problem on the square 
lattice, through a cell-to-cell estimate. For a fixed probability occupation p, we compute 
the specified wrapping probability as a function of the lattice size. Fig. HI shows R^^^ as a 
function of L. For p < pc, R^l^ decreases with increasing lattice size. For p > p^, it increases. 
Exactly at p^ one expects the wrapping probability to be independent of system size. Thus, 
by varying p in the critical region and searching for a constant R^l^ we obtain an estimate for 
Pc- The vertical dashed lines in Fig. E^a) represent the limits on the p values studied. From 
the data for R^^^ in Fig. H we have then p^ = 0.5928(2) and R^^^ = 0.523(4). Combining 
this estimate with the ones coming from the other quantities we obtain the values listed in 




FIG. 3: Maximum derivative of wrapping probabilities (-R^ for site percolation on the square 
lattice and i?^^ for site percolation on the hexagonal lattice) as a function of system size. Error 
bars are smaller than the symbol sizes. Linear regression of the data gives v = 1.334(4) for the 
square lattice and v = 0.866(1) for the hexagonal lattice. 
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FIG. 4: Wrapping probability as a function of the lattice size L for site percolation on square 
lattices. Different lines correspond to different concentrations p. The best constant horizontal line 
is represented by full triangles yielding = 0.5928(2) and i?^^ = 0.523(4). Error bars are smaller 
than the symbol sizes. 

Table im The results are quite close to the exact or expected ones, despite the small systems 
(sizes up to 128 x 128) and short runs (see table H)). Table ITTl also gives the results obtained 
from the present procedure to the bond percolation problem (which is easily implemented 
in the algorithm) with an excellent estimate of the known critical concentration. Moreover, 



the wrapping probabilities, within the error bars, are the same for site and bond problems, 
confirming the universal character of these quantities. 

Having demonstrated the good performance of the method in cases where exact results 
are available, we studied some three-dimensional lattices where data are not so ubiquitous 
in the literature. In particular, we treat the simple cubic, simple hexagonal and the HCP 
(hexagonal close-packed) lattices. To our knowledge, there are no results available for the 
wrapping probabilities on such geometries as well as no indication of their universal aspect 
regarding site and bond dilution. 

As an example, we show in Fig. ^c) and^d) the wrapping probabilities R^^^^ and R^^^ as 
a function of p for various lattice sizes, for site percolation on the simple hexagonal lattice. 
The corresponding scaling behavior of the derivative of R^^^ is depicted in Fig. |H1 an estimate 
for the critical exponent u may be extracted from these data. In Fig. Elwe show the estimate 

(x) 

for pc as well as the value of the wrapping probability Ri ■ The combined results are listed 



in Table IIIII together with the values for the simple cubic lattice and those obtained for 
the bond percolation on the HCP and simple cubic lattices. One can clearly see that the 
wrapping probabilities are distinct for different geometries, as is the critical concentration. 
Not only are our pc estimate comparable to the values obtained previously, but the critical 
exponents found here are close to the expected result for this universality class, namelv 

u = 0.8765(16)0 and z/ = 0.893(40)0 



u = 0.83(5) from series j2a], = 0.88(l)j2l, z/ = 0.8765(16)|2q and u = 0.893(40) j22l from 
Monte Carlo simulations. 

The data of Table IIIII are, up to our knowledge, quite new for these three dimensional 
lattices. A by-product of the present results concerns the universality of Rl{Pc) at the 
percolation threshold. One can clearly see that for site or bond percolation the wrapping 
probabilities of the simple cubic lattice are, within the error bars, the same (as well as for 
the problem in two dimensions depicted in Table ITT|) . This is in quite good agreement with 
the expected universal aspect previously obtained for the spanning probability in general 
dimensions and with both free and periodic boundary conditions 29|. 

Another interesting aspect of the three-dimensional lattices is the behavior of the quantity 
R^l\p) giving the probability of wrapping around two directions and not around the third 
direction. Fig. IHl shows such behavior for the simple cubic site diluted problem. In this case 
one can obtain an estimate of the critical concentration not only from the position of its 
maximum but also from the crossings at p^.. It is noted, however, that a more accurate value 
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FIG. 5: Wrapping probability as a function of the lattice size L for site percolation on 
hexagonal lattices. Different lines correspond to different concentrations p. The best horizontal 
line is represented by full triangles, yielding pc = 0.2624(4) and i?^^ = 0.429(5). Error bars are 
smaller than the symbol sizes. 

is achieved from the analysis of the crossings (as done in Fig. than from the location of 
its maximum. The same behavior occurs for the other lattice structures. 
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FIG. 6: Wrapping probability (p) as a function of p for the site percolation on the simple cubic 
lattice. From the position of the peaks one estimates pc = 0.3171(3) and from the procedure of the 
crossing region one gets pc = 0.3116(3). 



IV. CONCLUDING REMARKS 



We have seen that the results obtained using the present method are in good agreement 
with the exact (when available) or expected ones. It is important to note that this pro- 
cedure has been implemented using relatively small systems and short Monte Carlo runs; 
better results could be achieved in larger-scale simulations. In addition, from the computed 
wrapping probabilities at pc for the square and simple cubic lattices one can also confirm 
their universal aspect regarding site and bond dilution. The same should hold, of course, 
for the other lattice geometries. 

A slightly modified version of the present procedure, not using the cell-to-cell estimate, 
can also be applied to problems on two-dimensional lattices. Instead of tuning p, one can 
tune one of the Rl in the critical region. For a given quantity Rl (say, for example R^i^) we 
fix it at a specified value R* on the vicinity of the critical point, and proceed analogously 
to Newman and Ziff's original approach. Observe that R* is in this case a first estimate 
for Roo{Pc)- One can then compute, for each L, the intercept between the function Rl{p), 
previously evaluated, with the fixed R*. Each of these intercepts gives an estimate for p^, 
which is expected, for two-dimensional systems, to scale as L'"^ with cj) = 11/4. Therefore, 
plotting p^ X for different values of R*, we can estimate the true -Roo as well as the 

percolation threshold by looking for the value of R* that yields the best straight line. Table 
im reports the critical values so obtained for the two-dimensional lattice. The results are 
quite good and comparable to those obtained in the previous section for both the wrapping 
probabilities and the critical concentration. However, the errors are in general considerably 
larger. In three dimensions this procedure can be implemented only if we know the exponent 
beforehand. 
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TABLE I: Lattice sizes and run lengths (MCS) used in this work. The smallest and largest figures 
correspond to the total number of sites or bonds. In parenthesis we have the corresponding lattice 
size L. 

Lattice Smahest Largest MCS (10^) 

Site percolation 
Square 1 024(32) 16 384(128) 4.0 
Hexagonal 4 096(16) 21 952(28) 1.0-10.0 
Cubic 4 096(16) 21 952(28) 1.0-2.0 

Bond Percolation 
Square 2 048(32) 32 768(128) 2.0-4.0 
HCP 768(4) 20 736(12) 5.0 
Cubic 12 288(16) 52 728(26) 0.1-0.5 

TABLE IL Results for site and bond percolation on the square lattice. Errors in parenthesis affect 
the last digits. For each case, the first row shows the results described in section III and the second 
row those from the modified approach briefly discussed in section IV (in the latter method the 
exact value for the exponent v) is used. 



Two dimensions 



Lattice 


R{h) r{v) R{e) j^{b) ^ 


Pc(this work) 


Vc 


site 


0.517(4) 0.523(4) 0.692(3) 0.347(4) 1.333(5) 


0.592 7(1) 


0.592 T 




0.521(9) 0.524(6) 0.695(7) 0.353(5) 


0.592 9(3) 


0.592 1" 


bond 


0.521(2) 0.518(3) 0.691(3) 0.351(2) 1.331(3) 


0.499 95(15) 


1/2 




0.517(11) 0.519(13) 0.684(16) 0.348(6) 


0.499 8(4) 


1/2 


Exact'' 


0.5211 0.5211 0.6905 3516 4/3 






'^Ref.JlS]. 



TABLE III: Results for site(s) percolation on hexagonal and simple cubic lattices and for bond(b) 
percolation on simple cubic and HCP lattices. Errors in parenthesis affect the last digits. 



Three dimensions 

Lattice ii^^') R^'^ i?^'^) v Pc(this work) pc 

Hexagonal(s) 0.429(5) 0.332(5) 0.183(4) 0.467(6) 0.120(3) 0.867(14) 0.262 5(2) 0.262 3(2)^^ 

Simple cubic(s) 0.254(5) 0.255(5) 0.254(5) 0.456(7) 0.078(3) 0.877(12) 0.311 5(3) 0.311 6063(9)^ 

Simple cubic(b) 0.265(6) 0.266(6) 0.265(6) 0.471(8) 0.084(4) 0.868(11) 0.249 0(2) 0.248 8126(5)'^ 



HCP(b) 0.331(6) 0.443(6) 0.093(3) 0.561( 


f^0.052(3) 0.848(33) 0.120 3(2) 0.119 9(2)-^ 


"Ref. 


25]. 


''Ref. 


m. 


'^Ref. 





